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Non-spherical emulsion droplets can be stabilized by densely packed colloidal particles adsorbed 
at their surface. In order to understand the microstructure of these surface packings, the ordering of 
hard spheres on ellipsoidal surfaces is determined through large scale computer simulations. Defects 
in the packing are shown generically to occur most often in regions of strong curvature; however, 
the relationship between defects and curvature is nontrivial, and the distribution of defects shows 
secondary maxima for ellipsoids of sufficiently high aspect ratio. As with packings on spherical 
surfaces, additional defects beyond those required by topology are observed as chains or “scars”. 
The transition point, however, is found to be softened by the anisotropic curvature which also 
partially orients the scars. A rich library of symmetric commensurate packings are identified for low 
particle number. We verify experimentally that ellipsoidal droplets of varying aspect ratio can be 
arrested by surface-adsorbed colloids. 


Significance statement 

Emulsions, combinations of immiscible fluids such as 
oil and water, comprise a variety of commercially relevant 
systems, from ice cream to cosmetics. Individual droplets 
in these emulsions can be stabilized by the presence of 
colloidal particles and nonspherical droplets can be pro¬ 
duced by a mechanism of arrested relaxation, where par¬ 
ticles adsorbed at the droplet interface become crowded 
and obstruct further evolution of the surface toward the 
spherical ground state. The particles tend to have a high 
degree of crystalline order, but the curvature of the sur¬ 
face frustrates this order and introduces defects. In this 
paper we study the defect structures that form in pack¬ 
ings of hard particles on ellipsoidal surfaces, providing 
insight into the microstructure and stability of arrested 
systems. 


I. INTRODUCTION 

Emulsions—mixtures of two immiscible fluids—are 
ubiquitous systems with many applications in the food, 
oil, and cosmetics industries. At the microscopic level, an 
emulsion consists of droplets of one fluid embedded in a 
host fluid; the droplets are held in an equilibrium spher¬ 
ical shape by the interfacial tension between the two flu¬ 
ids. Emulsions with anisotropic droplets are of interest 
because for some applications, e.g. particle filtering in 
porous media[T], performance is improved with increas¬ 
ing aspect ratio. Anisotropic particles are also known 
to be more easily absorbed by cells, thus being effective 
as drug delivery systems [U [3]. Additionally, ellipsoids 
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Figure I: Sample packing oi N — 800 particles on a prolate 
ellipsoid of aspect ratio 2.6. (A) Side and (B) end views are 
shown; corresponding plots are shown for an oblate ellipsoid 
of aspect ratio 2.6 (C) from the top and (D) around the rim. 
Particles are colored by coordination number as computed 
from the Delaunay triangulation of the centroids — 5: brown; 
6: white; 7: blue; 8: light blue. 


fill space more efficiently than spheres [4], and through 
chemical functionalization, are a valuable component in 
the nano-architecture of hierarchical structures [5]. 

A mechanism for sculpting stable shaped droplets 
exists in Pickering emulsions, where the constituent 
droplets are stabilized by colloidal particles adsorbed at 
the interface|6]. The particles are strongly bound to 
the surface because they reduce the interfacial tension 
between the two immiscible phases[7]. Non-spherical 
shapes can be produced by a sequence of deformation, 
adsorption, relaxation and arrest as follows: Eollowing 
an initial deformation, for example by an applied electric 
field[8] or by the coalescence of two droplets[9]; during 
this process additional particles may become adsorbed 
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on the interface from the host fluid. The droplet then 
relaxes towards the equilibrium spherical shape, reduc¬ 
ing the surface area and causing the particles to become 
more densely packed. If the surface coverage of colloids is 
sufficiently high, they will become crowded and arrest the 
shape evolution of the droplet before a spherical shape is 
reached [9l [To]. 

The purpose of this paper is to identify the role that 
the anisotropic curvature present in an ellipsoid plays on 
the ordering of the particles. We assume the particles 
interact purely through volume exclusion. The quality of 
the packing of the final state, measured globally by cov¬ 
erage fraction as well as locally by coordination number, 
depends on this ratio: As T^jT^ 0, the particles are 
unable to rearrange themselves significantly and may get 
trapped in a glassy state, while for ^ oc, the relax¬ 

ation proceeds slowly and the situation resembles a clas¬ 
sical sphere packing problem. It is this latter quasi-static 
limit of the relaxation process that we shall examine in 
this work. 

Since the colloids are confined to a 2D surface, the 
arrested states tend to be quite crystalline as has been 
shown for spherical droplets or colloidosomes^^. These 
structures should, therefore, exhibit properties similar to 
2D elastic crystalline membranes [T2]-[22]. The presence of 
curvature frustrates the crystalline order and induces de¬ 
fects: particles which have more or fewer than six neigh¬ 
bors, and whose deviation from six-fold order can be 
quantified as a topological charge: particles with coor¬ 
dination number lower than six have positive charge and 
vice versa. Lone defects of positive or negative charge 
are known as disclinations. The topology of the droplet 
surface will determine the net defect charge, which is 12 
for a spherical topology [23]. Furthermore, there is a cou¬ 
pling of defects to the Gaussian curvature K. Because 
droplets with non-spherical geometries possess a varia¬ 
tion in Gaussian curvature along their surface, the de¬ 
fects should be non-uniformly distributed as theoretical 
studies have predicted [T8ll2Q] . 

In addition to the minimal number of defects required 
by topology, pairs of positive and negative defects called 
dislocations can occur. Droplets with a large system size 
, i.e. where the ratio R/r of the droplet size R to the 
particle size r is large enough, exhibit chains of defects 
known as scars pTl [T6] . For spherical droplets, a transi¬ 
tion has been shown: R/r is below a critical value only 
isolated defects occur. Above this ratio, scars appear and 
increase in length with R/r\W\. 

For surfaces of nonuniform curvature, the placement 
of the defects is an interesting question. The theory of 
curved elastic crystalline membranes [T4| predicts that de¬ 
fects and Gaussian curvature act as source terms in a 
biharmonic equation, 

V^x(f) = p(x) - K(f), (1) 

where x is a stress function and p is the defect charge 
density (a sum of point charges). The energy of such a 


system is, 

u = [ dAx{x){p{x) - K{x)), (2) 

Js 

which must be minimized with respect to defect number 
and defect position, with total defect charge conserved 
according to the surface topology. While this suggests 
that defect charges will be attracted to areas of like- 
signed curvature in order to minimize the source term, 
the fact that these systems are governed by a biharmonic 
equation suggests that the coupling of defects to curva¬ 
ture is nontrivial. This is in contrast to simpler ana¬ 
logues, for example electrostatics, governed by a Poisson 
equation. 

There are two important differences between an elastic 
crystalline membrane and a 2D arrested hard sphere sys¬ 
tem. First, in the hard sphere limit, the in-plane elastic 
constants of a hard sphere system are infinite. Second, 
arrested hard sphere systems are not able to explore their 
full phase space, and so one does not expect to find them 
in the optimally packed ground state. For these reasons, 
it should not be expected that a 2D hard sphere sys¬ 
tem is exactly described by the theory of 2D crystalline 
membranes, but due to the highly crystalline order that 
is exhibited, the behavior should be qualitatively simi¬ 
lar. Additionally, the relative wetting properties of the 
two fluids may induce a contact angle, leading to inter¬ 
particle interactions that may modify the ordering|24). 

In other systems in the packing limit, e.g. viral 
capsids [25] and small clusters of colloids [26], configura¬ 
tions with a high degree of symmetry are typically ob¬ 
served for certain special numbers of particles. Experi¬ 
mentally, these tend to be stable, and so the identification 
of possible symmetric packings may serve as a guide to¬ 
wards stable self-assembled micro-structures. We there¬ 
fore examine the packings systematically by aspect ratio 
a and particle number N to identify the symmetric con¬ 
figurations. 

In order to explore the role of surface anisotropy on 
the ordering of packed particles, we present the results of 
simulations of hard spheres packed onto ellipsoidal sur¬ 
faces using an inflation algorithm. Sample results are 
shown in fig. We investigate the effect of aspect ra¬ 
tio and particle number on the average distribution of 
defects on our surfaces and the structure of the defects 
themselves. We also identify highly symmetric configu¬ 
rations. Experimentally, we demonstrate that ellipsoidal 
droplets can be stabilized by surface-adsorbed colloids, 
and we compare the spatial distribution of defects in the 
experiments and simulations. Details of the model and 
simulations are presented in Methods. 

II. RESULTS AND DISCUSSION 

We employ an inflation packing algorithm in order to 
generate packings of spheres on ellipsoidal surfaces. The 
centroids of N equal sized spheres are bound to a fixed 
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ellipsoidal surface, either prolate or oblate, of aspect ratio 
a. The particles have hard-sphere interactions and diffuse 
as the particle radius is slowly incremented, until further 
inflation is precluded. Further details of the algorithm 
are given in Methods. 

Two sets of data were generated from which we ob¬ 
tained our results. One data set was used for studying 
the curvature-defect coupling and scar length, which con¬ 
sisted of packings with aspect ratio varying from 1.2 to 
4.0 in increments of 0.2 (for both the prolate and oblate 
cases: we consider the aspect ratio to be the ratio of 
the semi-major to semi-minor axis.) The particle num¬ 
ber was varied from 10 to 800 in increments of 10. Ad¬ 
ditional prolate packings were generated to study scar 
orientation, from aspect ratio 4.2 to 8.0 in increments 
of 0.2, from particle number 710 to 800 in increments 
of 10. 50 configurations were generated for each pair of 
parameters. The second data set was used for study¬ 
ing symmetry, where we are interested in lower particle 
numbers and a more fine-grained search of the parameter 
space. This data set consisted of packings with aspect 
ratio varying from 1.1 to 4.0 in increments of 0.1, and 
particle number varying from 3 to 200 in increments of 
unity. 80 configurations were generated for each pair of 
parameters. 


A. Defect Distribution 
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Figure 2: Defect number density for (A) prolate and (B) 
oblate ellipsoids of varying aspect ratio: blue is 1.2; yellow 2.6; 
purple 4.0. Note the small secondary peak near zj zq — 0.4 at 
a = 4 in the prolate case. Example configurations of a = 4 are 
shown as insets. Defect charge density is shown for (C) pro¬ 
late and (D) oblate ellipsoids of a = 4. The green points rep¬ 
resent the net charge density, and the brown and blue points 
represent the density of positive and negative defects, respec¬ 
tively. The secondary peak in (A) is also visible in the positive 
and negative charge densities in (C). In (D), there is a net 
negative defect charge density near z/zo = 0.4, despite the 
Gaussian curvature being positive. Lines are guides to the 
eye. 


We first examined the distribution of the defects as a 
function of the aspect ratio. Defect locations were deter¬ 
mined by assigning a defect charge q = 6 —c to each par¬ 
ticle, where c was the coordination number determined 
from the Delaunay triangulation of the particle positions 
(see Methods). The surface was partitioned into equal- 
area axisymmetric regions and the number of defects in 
each region counted. Each segment has a different aver¬ 
age Gaussian curvature with regions near the poles hav¬ 
ing larger curvature for prolate and the reverse for oblate 
ellipsoids. In fig. for prolates and fig. for oblates, 
the defect number density is shown as a function of the 
axial position zj averaged over the ensemble of simula¬ 
tions at fixed aspect ratio and particle numbers ranging 
from 710 < A/" < 800. Generically, it is apparent that 
defect number density increases with the Gaussian cur¬ 
vature, as expected. For prolate ellipsoids at low aspect 
ratio, the defect number density increases monotonically 
with respect to K. At higher aspect ratios, there is a 
small secondary peak in segments with low Gaussian cur¬ 
vature. We verified this occurs for other ranges of particle 
numbers N > 210. 

In order to understand this, we plot separate defect 
charge densities for positive and negative defects in fig. 
[2](7, as well as the net defect charge density. The anoma¬ 
lous peak is apparent in both the separate positive and 
negative defect charge densities, but not in the net de¬ 
fect charge density, indicating that the excess defects are 
taking the form of neutral dislocations or scars. 


In fig. 1^, we see that for oblate ellipsoids, the de¬ 
fect density again increases near the more highly curved 
regions. Fig. W reveals, however, that the coupling be¬ 
tween defect charge and curvature is again complicated: 
while there is a peak in positive defects at the highly pos¬ 
itively curved edge of the surface, there is a high density 
of negative defects surrounding this, and the net defect 
charge density is actually negative near z/zq = 0.4. 

These results display a non-trivial interaction between 
defects and curvature. While the regions of highest Gaus¬ 
sian curvature contain the highest density of defects, 
the defect density is not a simple monotonic function of 
Gaussian curvature. This is apparent in the defect num¬ 
ber density in the prolate case, and in the defect charge 
density in the oblate case. The fact that the defect charge 
density can be negative in regions of positive Gaussian 
curvature is especially surprising. However, this is not 
necessarily inconsistent with eqs. [^andj^ which imply 
complex defect behavior. Further investigation is war¬ 
ranted to confirm whether the continuum elastic theory 
gives results similar to the hard sphere packings here. 


B. Scar Orientation 

We next determined whether the scars are oriented by 
the curvature anisotropy of the surface. To do so, we 
consider a local scar orientational distribution function 
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Figure 3: Orientation of the scars relative to the curvature 
anisotropy. (^) A configuration with a typical scar. (B) 
Close-up of the scar. Black lines show edges in a graph com¬ 
prising the scar. The red dashed line shows a chain of length 
3. Results are shown for {C-E) prolate and (F-H) oblate 
ellipsoids. The C 2 {D, G) and C 4 {E, H) order parameters 
for prolate and oblate ellipsoids, respectively, are plotted as 
a function of aspect ratio for different regions along the sym¬ 
metry axis of the ellipsoid: green corresponds to the center, 
orange to the mid-region, and blue to the ends. (C) and (F) 
show the ODF of chains in the center, mid-regions, and ends 
of the ellipsoid, respectively, for prolate ellipsoids of aspect 
ratio 8 in (C) and oblate ellipsoids of aspect ratio 4 in (F). 
Insets of (C) and (F) illustrate the regions used for spatial 
binning. 


(ODF) f{a) where the angle a is measured locally in the 
tangent plane relative to the uniaxial axis of the ellipsoid. 
The ODF may be expanded as a Fourier series, 

fia) = E CnCos{na). (3) 

n 


The average value of the first two non-zero coefficients, 
C 2 = {cos{2a)) and C 4 = (cos(4a)), were calculated 
for our ensemble of packings. These quantities are or¬ 
der parameters for orientational order as they vanish if 
the scars align isotropically with the curvature. C 2 quan¬ 
tifies nematic order, i.e. uniaxial orientational order and 
C 4 quantifies quadrupolar order. 

To determine the scar orientation, we studied contigu¬ 
ous chains of defects as shown in fig. and B. Given a 
packing and its Delaunay triangulation, the neighboring 
defects around each defect are identified. These adja¬ 
cent pairs become the edges of graphs of contiguous de¬ 
fects. Two defects are identified as the ends of a chain of 
length I if they are within a connected graph of defects 
and the shortest path between them contains I edges. 
Once a chain of length I is identified, its orientation rel¬ 
ative to the local principal directions—i.e. the polar and 
azimuthal tangent vectors te and tcj), respectively (see 
eq. 0 in the Appendix for the parametrization of the 
surface)—is calculated thus: given a pair of chain end¬ 
points, their separation vector is projected onto the sur¬ 
face at each endpoint, giving components along tg and 
These components are then averaged between the end¬ 
points, and the angle a that the resulting vector makes 
with to is recorded as the orientation of the chain. The 
^-component of the midpoint of each chain is recorded 
as its position and is used to examine how the coupling 
varies across the surface. 

The analysis was applied to an ensemble of simulation 
results as follows: For a given aspect ratio, the orienta¬ 
tions of all chains of length I are collected across simula¬ 
tions with N G [710,800] in increments of AA^=10 (with 
50 results at each N resulting in 500 simulations). Or¬ 
der parameters C 2 and C 4 are then calculated from this 
ensemble. Because the curvature anisotropy varies with 
the 2 :-coordinate along the surface, results can be divided 
according to their position. In our analysis, we exclude 
scars in the regions near the poles which make up 10 % of 
the surface area as here the curvature tensor is degenerate 
and the alignment is undefined. The rest of the surface 
is broken into six equal-area, azimuthally symmetric re¬ 
gions, as illustrated in the insets of fig. C and F, and 
data from symmetric regions on opposite halves of the 
ellipsoid are combined. A chain length of 1=3 was used 
as this is long enough to capture scar behavior while hav¬ 
ing enough chains for statistical purposes. Shorter chain 
lengths show a weaker tendency to orient. 

The behavior exhibited by prolate ellipsoids is rather 
complicated, as seen in the plots of order parameter ver¬ 
sus aspect ratio in fig. D and E. In the center region 
near the equator, scars are nematic along the to direction 
between aspect ratio 3.6 and 6 . At higher aspect ratio 
this center region is very flat, leading to fewer scars, and 
so any orientational order is insignificant. In the mid¬ 
regions between the equator and poles, scars become ne¬ 
matic along the direction at aspect ratio 4.4, and then 
transition to nematic along the to direction at aspect ra¬ 
tio 6.4. Scars near the poles show nematic order along 

















5 


to above aspect ratio 2, although this order peaks near 
aspect ratio 5, then drops to C2=0 at aspect ratio 6.4 
before increasing again. Interestingly, scars on highly 
prolate ellipsoids can also show C 4 , order. This appears 
in the mid regions above aspect ratio 5.2, and in the end 
regions above aspect ratio 6. 

The chain ODFs for prolate ellipsoids of aspect ratio 
a = 8 in fig. HP illustrate the trends that appear at 
high aspect ratio. It is apparent from the green curve 
that that there are few chains in the relatively flat center 
of the ellipsoid. The orange curve shows a high degree 
of nematic order directed along the polar direction in 
the mid-region, and the blue curve for the ends shows 
nematic order along the polar direction, as well as a peak 
between the directions of principal curvature, which is 
indicative of negative order. 

The case of scar orientation on oblate ellipsoids is more 
straightforward. The order parameters are plotted as a 
function of aspect ratio for different azimuthally symmet¬ 
ric regions across the surface, in fig. G and H. Scars at 
the equator exhibit a high degree of nematic order in the 
t(f) direction, which increases linearly with aspect ratio 
up to a = 4. This is unsurprising, because the curvature 
on highly oblate ellipsoids is localized to a nearly one¬ 
dimensional region around the equator of the ellipsoid, 
and so one expects the scars to form there, aligned along 
the equator. There is also a small degree of C 4 ordering. 
In the regions midway between the equator and poles, 
there is a weak coupling of scars along the to direction. 
These trends are illustrated for a = 4 in fig. HP :the green 
curve for the edges displays a peak near the azimuthal 
direction, whereas the orange and blue curves show that 
there are fewer chains without much order in the flatter 
regions. 

While the scar orientation results for the oblate case 
are easily understood, the ordering of the scar orientation 
on prolate ellipsoids is far more complicated. The ori¬ 
entation varies greatly depending on chain position and 
ellipsoid aspect ratio. Especially surprising is the emer¬ 
gence of C 4 ordering, which corresponds to a tendency 
for chains to align in a direction intermediate to the di¬ 
rections of principal curvatures. 


C. Scar Transition 

As is well known from previous work |14l I16| . packings 
of spheres on spherical surfaces exhibit a transition: For 
low particle numbers, only the twelve defects required by 
topology are present; above a critical particle number Nc^ 
it is favorable for larger defect structures to occur, typ¬ 
ically chains of scars extending from a core disclination. 
Increasing N above Nc leads to a monotonic increase in 
average scar length. 

From our simulation results of packings with 10 < A" < 
800, we calculated the average number of excess dislo¬ 
cations per topologically required disclination for each 
(a. A). Defects were weighted in the analysis by the ab¬ 



Figure 4: The number of excess dislocation defects per scar 
on (A) prolate ellipsoids and (B) oblate ellipsoids. For low 
aspect ratio near 1, there is a clear scar transition, which 
is not present at aspect ratios far from 1. The inset in (B) 
shows a highly commensurate oblate packing with A = 140 
and a = 2.6. Note that data for oblate ellipsoids with A = 10, 
a > 2.0 and A = 20, a > 3.0 has been excluded. . 


solute value of their charge. Given that there are two 
disclinations per dislocation, and 12 core disclinations, 
the number of excess dislocations per scar is calculated 
thus. 




( 4 ) 


where the sum is taken over all defects. This quantity 
captures the same information as the scar length but is 
easier to calculate, as individual scars are often not well 
defined. 

Results of the analysis are displayed in fig. Prolate 
ellipsoids [fig. [^] show the experimentally observed be¬ 
havior for low aspect ratio: for A < 100 particles there 
are few excess defects, but at higher particle numbers 
there is a roughly linear increase in the number of excess 
defects. As aspect ratio increases, however, the transi¬ 
tion is softened such that there is a smooth increase in 
excess defects with A. This is reminiscent of how applied 
fields soften phase transitions [57]; here the anisotropy of 
the curvature seems to play a similar role. 

The oblate packings show the same trends [fig. 

There is, however, an additional feature that stands out. 
At A = 140, a > 2, there is a set of nearly scar-free 
configurations. This is due to commensurability as the 
particle number and surface geometry for these cases are 
compatible with a highly symmetric packing with only 
the minimally required defects, as seen in the inset of 
fig. HlB. Similar commensurability issues occur in other 
systems, e.g. sphere packings on cylinders [28]. 

A striking difference between these results and those 
from a previous study is that here, for hard particles, 
the transition occurs at a lower particle number; in ref. 
m it was seen at Ac ~ 400 using colloidal particles 
with a soft repulsive interaction. We therefore performed 
simulations (see Methods) using two different potentials, 
V = d~^ and V = d~^ (where d is the interparticle sepa¬ 
ration), the results of which are shown in fig. For soft 
particle packings, we take the average scar length of the 
five lowest energy configurations obtained out of an en¬ 
semble of 50. For the hard spheres, Ac ~ 80, while for the 
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Figure 5: Excess dislocations per scar as a function of particle 
number for hard (blue) and soft V — Xjd (orange) and V — 
1/d® (red) interactions. Inset (A) is a hard particle packing 
and inset (B) is a soft particle packing. The arrow indicates 
the particle number of the inset packings. 


two soft potentials the transition occurs around Nc ~ 200 
(which appears to be within the uncertainty of the result 
presented in ref. M)- The defect number increases at 
the same rate with respect to particle number for both 
soft potentials. This supports the conclusion in ref. m 
that, for soft particles, the scar transition does not de¬ 
pend on the specific form of the particle potential. For 
hard particles we have quantitatively different behavior. 
Visual inspection of hard and soft sphere configurations 
reveals that hard sphere configurations possess gaps (fig. 
[^) . It is rare to find a lone disclination; it is much more 
common to find a disclination attached to one dislocation 
(i.e. a small 5-7-5 scar) adjacent to a gap in the packing. 
This isn’t seen in soft particle configurations (fig. iB), 
as the energy penalty is too high, rather a particle can 
be squeezed to fill in the gaps. The fact that hard par¬ 
ticle packings tend to have gaps makes them especially 
suitable for chemical functionalization as described in ref. 

m 


D. Packing Fraction and Symmetry 


We now turn to how the packing fraction varies with 
respect to particle number and ellipsoid aspect ratio. 
To simplify the calculation we make the approximation, 
valid for large V, that the area covered by a particle is 
its projection onto a flat 2D surface. 


(j) = 


N7Tr‘^ 

A 


( 5 ) 


where A is the area of the underlying surface. We 
checked the validity of this estimate by numerically inte¬ 
grating the area of intersection between the surface and 


the spheres on oblate surfaces of aspect ratio 4.0, and 
found that the difference between our estimate and the 
true value is very small: using the projected area un¬ 
derestimates the packing fraction by approximately 1% 
for packings with V = 100 and 0.1% for packings with 
N = 800. 

For large V, the packing fraction increases slightly 
with aspect ratio. This is because for large a the cur¬ 
vature—and hence the defects—are mainly localized to 
the poles on prolate surfaces or the equator on oblate sur¬ 
faces and so more of the surface can be covered by the 
planar hexagonal packing, consistent with the results of 
the above subsections on the Defect Distribution and Scar 
Transition. For low V, the opposite tends to be true; the 
packing fraction decreases with aspect ratio. However, 
the trend is more complex and the packing fraction is 
sensitive to both N and a at low N. Visual inspection 
of these configurations reveals that for specific combina¬ 
tions of N and a, the packings have a high degree of 
symmetry, suggesting a commensurability effect, such as 
that seen in the Scar Transition subsection above. 

To identify these commensurate combinations, we con¬ 
ducted a more thorough search for symmetric packings 
using the second data set. An arbitrary packing must 
break the ellipsoidal symmetry group of the surface and 
hence must belong to some finite subgroup of D^^] most 
packings at high particle number do so trivially, retaining 
only the identity element. Defining a suitable inner prod¬ 
uct (A, B) that measures the distance between two pack¬ 
ings, a packing possesses a symmetry C if {A^CA) = 0 
where C is a group element of Dooh- The elements C 
can be constructed from the group generators: i) an in¬ 
finitesimal rotation about the ellipsoid symmetry axis; ii) 
spatial inversion, and iii) a rotation by tt about an axis 
perpendicular to the symmetry axis. 

We used a norm (A, B) defined such that. 


(AS) 



(6) 


where the di and hj are the positions of particles in pack¬ 
ings A and B, respectively: for each particle in A, the 
closest particle in B is found and the separations be¬ 
tween these pairs are divided by the particle radius. The 
root mean square of these normalized separations is then 
taken as the inner product. From this, together with the 
group generators, all symmetries such that (A, CA) < e,a 
threshold separation were found. From this catalog of 
symmetries, for a particular configuration the appropri¬ 
ate group was determined. From a collection of config¬ 
urations with a given (V, a), the most symmetric con¬ 
figuration was chosen by the following procedure. First, 
the configurations with the largest symmetry group were 
identified. Then, for each of these configurations, the 
symmetry group element with the highest symmetry 
norm was identified. Finally,the configuration with the 
minimum highest symmetry norm was chosen as the most 
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Figure 6: The symmetry landscape for packings with vary¬ 
ing particle number and aspect ratio, using a symmetry norm 
cutoff of 0.1. (A) shows the chirality and the order of the 
largest symmetry group found. Orange represents chiral pack¬ 
ings and blue represents achiral packings. The boldness of 
the color corresponds to the order of the packing’s symmetry 
group as shown in the key. Note that packings whose only 
symmetry is the identity are colored white to distinguish them 
as being trivially symmetric. Sample packings are shown: B) 
a chiral packing with N = 74, a = 2.5; C) an achiral packing 
with N = 74, a = 1.5 — note that (B) and (C) have the same 
particle number, but show different chirality for different as¬ 
pect ratio; D) a packing with fourfold rotational symmetry 
with ^" = 69, a=1.4;E)a packing with fivefold symme¬ 
try N = 76, a = 2.4. Light brown particles have c = 4. (F) 
shows the degree of rotational symmetry of each configuration 
about its ellipsoidal symmetry axis. Note that for both (A) 
and (F), no data is shown for a = 1 (spheres) as the spherical 
symmetry group is not a subgroup of Dooh- Sample packings 
are shown for G) N = 30, a = 2.4; F) iV = 34, a = 2.5; I) 
N = 38, a = 2.7; J) N = 46, ; these packings all occur in the 
diagonal band of fourfold rotational symmetry in the top left 
of (F). 



symmetric. 

The results of this analysis are displayed in fig. |§4 
showing the order and chirality of the symmetry group 
of the best packing for each combination of particle num¬ 
ber and aspect ratio. The degree of rotational symmetry 
for each packing is shown in fig. HF- One striking feature 
is that, for certain particle numbers, long vertical stripes 
appear in the plots representing commensurate aspect ra¬ 
tios for that particle number. Furthermore, low N favors 
achiral packing while chiral packings occur more often for 
higher particle number. For prolates the stripes occupy 
a narrow range of aspect ratio and occur in band-like se¬ 
quences described by a straight line a = mN with slope 
m. Each of these sequences corresponds to a different de¬ 
gree of rotational symmetry and the particle numbers 
in the sequence are separated by Inspecting the con¬ 
figurations in a single sequence, the difference between a 
configuration with N particles and the next with N -\-nr 
particles is that an additional row of particles has been 
inserted in the space created by the longer aspect ratio. 
This is illustrated by a sequence of configurations with 
fourfold rotational symmetry in fig. G-J. 

For oblate ellipsoids, the symmetric configurations for 
N particles occur at a much broader range of aspect ra¬ 
tios and symmetric configurations are observed at much 
higher N and tend to have six-fold rotational symmetry. 
The reason for this is that the high curvature at the end 
of the prolate ellipsoids accommodates n^-fold defects at 
the poles, and these appear to determine the rotational 
symmetry for the entire configuration; for oblates, the 
poles have low curvature and promote hexagonal pack¬ 
ing, hence causing six-fold rotational symmetry to be 
more common. Interestingly, other degrees are present 
including = 4 and = 5 and these configurations 
contain regions of highly oblique packings (fig. D and 
E). 

In general, these symmetric packings are notable be¬ 
cause they contain a high degree of hexagonal ordering 
over much of their surface, with evenly spaced defects 
throughout. This high degree of regularity should pro¬ 
vide stability to the packed structure, and reduce the 
likelihood of failure from irregularly spaced defects. 


III. EXPERIMENT 

An experimental realization of ellipsoidal arrested 
droplets was performed to confirm the stability of these 
structures. Ellipsoidal droplets with arrested interfaces 
are produced by preparing a Pickering emulsion and then 
mixing the emulsion to deform and arrest the droplets in 
an elongated shape. Details are given in Methods. 

Fig.0 shows several examples of the arrested droplets 
observed. Because the curvature of a droplet is signif¬ 
icant across its surface, several focal planes have been 
combined in the images in order to study the packing of 
spheres on the drop surface. Particle coordinates are de¬ 
termined by finding the local brightness maxima in the 














Figure 7: Experimental data for particle-stabilized droplets 
of aspect ratio (A, B) 1.6, (C, D) 5.1, and (E, E) 3.0. Scale 
bars represent 15//m. (A, (7, E) Microscope images; (E, E, 
E) Reconstructed particle positions, colored by coordination 
number as determined by Delaunay triangulation of the par¬ 
ticle centroids — 4: light brown, 5: dark brown, 6: white, 
7: dark blue, 8: light blue, 9: purple. In general, defects 
are more common and are more likely to be found at low- 
curvature regions of the droplet in the experiments than in 
simulations. 


image, recording their coordinates, and correcting for any 
unrealistic results via direct comparison with the exper¬ 
imental images. 

Arrest is able to preserve shapes identical to interme¬ 
diate states of droplets in an elongation field [29], as seen 
in the fig[^ and[7](7, and even shapes resembling sections 
of such shapes as in the case While the dynamical 

formation of these shapes was not studied, it is clear that 
a wide range of geometries can be formed. We note that 
the droplet of aspect ratio 5.1 has a spherocylindrical 
geometry, as opposed to ellipsoidal. 

Fig. 0 E, E, and F shows the results of a Delaunay 
triangulation of the sphere coordinates. We do not dis¬ 
play particles at the boundary of the triangulation, as 
they include spurious edge defects identified as a result 
of the boundary rather than the ordering of the parti¬ 
cles. In each case the arrested state of the interfacially 
adsorbed spheres is evident from the visible regions of 
crystalline order. Generally, however, the experimental 
droplets contain more defects than the simulated pack¬ 
ings. In fig. a high degree of hexagonal close-packing 
is noted near the ends of the droplet, while the center 
of the structure is more disordered with a higher defect 
density. Three important factors present in the experi¬ 
ment that are not accounted for in the simulation may 
contribute to this. First, the evolution of the surface as 
it relaxes will influence particle rearrangement. Differ¬ 
ent parts of the surface will grow or shrink at varying 
rates, affecting where crowding first occurs. Second, par¬ 
ticles adsorbed at an interface will not act as purely hard 
spheres. Capillary interactions caused by the deforma¬ 
tion of the surface by the particles will lead to attractive 


interactions between particles[21|. This may lead to ag¬ 
gregation of particles during relaxation and is likely to 
influence the final ordering of the arrested state. Finally, 
as discussed in the introduction, the experimental relax¬ 
ation does not take place quasistatically, as is posited 
by studying the packing limit; it is highly likely that 
the particles are arrested in a nonoptimal and possibly 
metastable glassy state. A study of the role of these dy¬ 
namical influences on the order is in preparation. 


IV. CONCLUSION 


In this paper, we show that defects in the packing of 
hard spheres onto an ellipsoidal surface couple nontriv- 
ially to the curvature. For low aspect ratios, the defects 
occur at regions of high curvature as predicted by pre¬ 
vious studies; additional secondary peaks in the defect 
distribution occur in less-curved regions for prolate el¬ 
lipsoids of sufficiently high aspect ratio. As previously 
observed for packings on a spherical surface, above a crit¬ 
ical particle number the defects take the form of chains 
or “scars” rather than isolated defects. This scar tran¬ 
sition occurs at a lower particle number than the previ¬ 
ously studied case for soft interparticle interactions, and 
is softened by the presence of anisotropic curvature. The 
alignment of the scars with the curvature is more com¬ 
plicated: in flat regions, there is no alignment; in inter¬ 
mediate regions, there is weak uniaxial alignment with 
the minimum curvature; in regions of strong curvature, 
quadrupolar alignment is seen. We identified a rich cat¬ 
alog of symmetric configurations from our simulations, 
each belonging to a subgroup of the ellipsoidal symme¬ 
try group. Plotting the subgroup order in (V, a) space 
reveals commensurate surfaces that promote symmetric 
packings. Finally, we were able to use the mechanism of 
arrest to sculpt ellipsoidal Pickering emulsion droplets of 
varying aspect ratio, demonstrating the validity of the 
fundamental idea. While careful analysis of these ex¬ 
perimental packings reveals scars as predicted, the de¬ 
fects appear to agglomerate in regions other than those 
of strongest curvature, suggesting that dynamical effects 
play a significant role in the ordering as well as the geo¬ 
metric effects studied here. 
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Methods 


of the conformal factor, 


Experimental preparation of arrested droplets 

Emulsions are first prepared by mixing 3% w/w 
monodisperse 1.5 pm diameter precipitated silica par¬ 
ticles (Nippon-Shokubai KE-P150) into hexadecane 
(Sigma-Aldrich, 99%) [9]. A volume of the silica- 
hexadecane dispersion is then emulsified into an equal 
volume of deionized water by manual shaking for three 
minutes. The emulsion was then aged for 24 hours 
and inspection revealed a small fraction of elongated 
droplets. Imaging of the droplets is carried out on a 
Leica DM2500M light microscope using phase contrast 
optics. 

Hard-sphere simulations 


u{0) = f yW)dO\ (9) 

J TV j2 

which can be inverted to find 0{u). We do an approxi¬ 
mate inversion by calculating u{0) for values of 0 from 
0 to TT in increments of tt/IOO and using a high order 
polynomial least squares fit on these points, enforcing 
equality between the fit and exact values at the endpoints 
^ = 0 and 0 = IT. The conformal coordinate v is simply 
v{4>) = <i>. 

Given the definitions above, diffusion steps are taken as 
follows. An unsealed step size is chosen for each direction, 
and Ar’o, from a normal distribution with variance 
1. These are scaled by the simulation step size a and by 
the inverse of the conformal factor to give step sizes in 
the (i4, v) conformal space: 


We employ a stochastic inflation packing algorithm in¬ 
spired by the Lubachevsky-Stillinger algorithm, which is 
known to yield packings of high coverage fraction [30]. In 
each packing simulation, a fixed ellipsoidal surface, either 
prolate or oblate, is chosen with aspect ratio a and the 
length of the semi-minor axis is fixed to be unity in di¬ 
mensionless units. Particles are modeled as monodisperse 
hard spheres of radius r that is slowly increased during 
the simulation. The number of particles N is specified 
and particles are deposited at the start of the simulation 
by random sequential adsorption such that the center of 
each particle is constrained to lie on the surface of the 
ellipsoid. Initially, r is such that the packing fraction is 
(j) = 0.05. 

The algorithm proceeds by two kinds of moves: i) 
Monte Carlo diffusion steps where particles are moved 
randomly along the surface and ii) inflation steps where 
the radius of all particles is increased by 6 r. In each 
diffusion step, N individual Monte Carlo moves of ran¬ 
domly chosen particles are attempted. The step size is 
chosen randomly using a Gaussian distribution, as de¬ 
scribed below. Only moves that do not result in overlap 
are accepted, with overlaps checked for in the 3D config¬ 
uration frame. 

The moves are performed in the 2D space of confor¬ 
mal surface parameters hence yielding a radially 

symmetric probability distribution of moving a certain 
arclength s in any tangential direction from the current 
location. The surface is parametrized as, 

x{0, (j)) = (xo sin 0 cos 0, xq sin 0 sin 0, zq cos ^) , (7) 

where xq = I, 2:0 = a for prolate surfaces and xq = a, 
Z{) = 1 for oblate surfaces. The determinant of the metric 
is. 


9 {d) = txo sin(6»)2 {z^ + a;o + (zq - xl) cos(26»)) , (8) 

and the conformal parameter u is given by the integral 


Au 


Av 


aAuQ 


V9{d{u)) 

aAvo 

V9{d{u)) 


( 10 ) 

( 11 ) 


These steps are used to update the previous u and v 
coordinates of the particle, which are then transformed 
to the 0 and (j) coordinates as explained above. Einally, 
the surface parametrization eq. [^is used give the particle 
coordinates in the 3D configuration space. 

Because 0 must have a value between 0 and tt, we take 
the following step if it falls outside this range at any 
point. If uis greater than u{0) (less than u{7t))^ we set 
u = 2u{0) — u {u = 2u{tt) — u) and v = mod ('i; + tt, 2tt) , 
i.e. we allow the particle to pass over the coordinate 
singularity at the poles of the surface. 

As the particles diffuse, a is varied in order to more 
efficiently explore relevant areas of configuration space 
(leading to large steps when the configuration is loosely 
packed and smaller, more relevant steps when tightly 
packed.) The initial value of a scales with the square root 

of the ellipsoid surface area A, cFinit = I x 10“^ Af¬ 
ter each time step, the fraction of attempted moves that 
were accepted is calculated. The length scale a is then 
decreased by 1% if the acceptance fraction is < 0.5 and 
increased by 1% otherwise; a is reset after each inflation 
(described below) to its initial value. Bounds are im¬ 
posed such that I X 10“^ < cr < I. Adjusting a leads to 
improved performance of the algorithm as the particles 
can diffuse more when they are less densely packed and 
take smaller steps (which are more likely to be accepted) 
when they are more densely packed. We do this as it 
is known that adaptive algorithms lead to packings of 
higher density [31]. We emphasize that in this work the 
Monte Carlo approach is used as an optimization strat¬ 
egy; it is not intended to, and indeed cannot, replicate 
the physical process by which the structures form since 
the cr updates are non-Markovian. 
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After M = 100 diffusion steps, an inflation step is 
performed where the particle radius is increased slightly 
(“inflated”) either by a specified fixed amount Ar = 

1 X 10“^ or by the half of the largest amount al¬ 
lowed that would not result in the overlap of any pair of 
particles, whichever is smaller. 

The halting criteria for these simulations is as follows: 
every L = 100 inflation steps, the relative change in cov¬ 
erage fraction A0 is calculated. If this is less than a 
specified value Acjyfoi = 10“^ then the simulation halts. 


Soft particle simulations 


In order to compare our results regarding scar forma¬ 
tion in hard particle packings to other work involving 
particles with soft interactions, we performed a set of sim¬ 
ulations using a modified Monte Carlo algorithm which 
incorporates a soft interparticle potential. In order to 
test potentials of different softness, the interparticle po¬ 
tentials are set as either Uint = d~^ or Uint = (where 
d is the center-to-center distance between particles). The 
particles diffuse similarly to the hard particle simulation 
with two differences: the average step size a is constant 
for all moves, and moves are accepted or rejected using 
a Metropolis scheme [32]. with acceptance probability 


P = 


exp(-AP/A:BT) 


AP < 0 
AP > 0 


( 12 ) 


where AU is the change in the system energy after a sin¬ 
gle particle move. The initial temperature is set by using 
a rough estimate of what the energy of a single particle 
in the final configuration will be assuming six-fold or¬ 
dering and that nearest neighbor interactions dominate: 
To = 6f/j„((2rest)/fcB, where Test = \/A/N is an esti- 
mate of the average particle separation. The system is 
annealed by multiplying the temperature by 0.99 after ev¬ 
ery 100 sets of diffusion moves until exp(—AP/Z cbT) ^ 1 
within machine precision. After every 100 sets of diffu¬ 
sion moves, the change in energy is recorded and the 
simulation halts once this change in energy is less than 
1 X 10-^^ 


Defect analysis of simulations 


To analyze defects in the simulated configurations, we 
use a ball-pivoting akorithm[33] in the mesh-generation 
software Meshlab to generate triangulations of the parti¬ 
cle centroids. These triangles are then equiangulated by 
a custom script to remove narrow triangles. Edges are 
flipped in random order and accepted if they improve the 
triangulation; this is repeated until a full sweep of the 
mesh yields no further improvements. From these op¬ 
timized triangulations, the coordination number of each 
particle is given by the number of particles to which it is 
connected. 
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